#ifndef S2DMESH_H
#define S2DMESH_H

#define sigma 100
#define dExpB 0.02
#define thetaExp 1.5707963267949
#define sorFac 0.9
#define resMin 1e-5
#define resMinPQ 1e-4
#define updateFac 0.5

#include "s2dMesh.h"

#include "mesh.h"
#include "readFile.h"
// #include "diff2d.h"
#include <vector>

using namespace std;

namespace s2dMesh{
    extern int xNum,yNum;
    extern vector<vector<mesh::point2d> > curveList;
    extern double boundaryRefineB;
    extern double boundaryLayer;
    static double **xMesh;
    static double **yMesh;
    int getS2dMeshPara(vector<fileData> meshPara);
    int getCurvePoint(vector<mesh::point2d> wallPts,int num,bool straight);
    int meshInit();
    int writeMesh();
    int writeMeshTemp();
    int freeMem();
}

namespace s2dMeshOpt{
    static double **aCoef;
    static double **bCoef;
    static double **cCoef;
    static double **qSource;
    static double **pSource;
    static double **xMesh;
    static double **yMesh;

    int optMeshInit();
    int optMesh();
    int writeResidual();
}

namespace s2dMeshOpt{
    static int iter=0;
    static double pqResInit=0;
    static double **xMeshNew;
    static double **yMeshNew;
    static double **phi;
    static double **psi;
    static double **xKsi,**xEta,**yKsi,**yEta;
    static double **xKsiKsi,**xEtaEta,**yKsiKsi,**yEtaEta;
    static double **j;
    static vector<double> residual;

    int s2dMeshGS();//done
    int copyBdryXY();//done
    int copyXY();//done
    int calcABC();
    int calcBdrABC();
    int calcItnABC();//done
    // int calcBdrABC();
    int calcBdrPQ();//to
    int calcPQ();
    int calcPQInit();
    // int calcBdrJ();
    int calcJ();
    int calcBdrPhiPsi();
    int calcPhiPsi();
    int calcBdrPhiPsi0();//done
    int calcBdrDiff();
    int calcItnDiff();
    int calcBdrDTheta();
    int updateBdrPoints();
}

namespace basic{
    vector<mesh::point2d> interpCurve2d(vector<mesh::point2d> curve,vector<double> x,bool interpByX);
    mesh::point2d interpPoint2d(vector<mesh::point2d> curve,double x,bool interpByX);
    vector<double> linspace(double min,double max,int num);
    vector<mesh::point2d> xyVecToCurve(vector<double> x,vector<double> y);
}

#endif